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ABSTRACT 

A method for tracking neutrons through a slab penetrated by a helical voided duct is 
presented and verified. FORTRAN listings of computer subroutines that perform the 
numerical analysis are included; these subroutines are designed to be used with FASTER, 
a high-efficiency Monte Carlo code . Neutron number spectra were calculated at several 
points in and near the emergent end of a duct through a 20- and a 100-centimeter thick 
water slab having a point source at the center of the inlet. 
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HELICAL DUCT GEOMETRY ROUTINE FOR THE SHIELDING 
COMPUTER PROGRAM ’’FASTER" 
by Thomas M. Jordan* and Millard L. Wohl 
Lewis Research Center 

SUMMARY 

A method for tracking neutrons through slabs penetrated by helical voided ducts is 
presented and verified by means of the FASTER computer code. A point fission source 
is located at the center of the entrant mouth of the duct. Neutron number spectra com- 
puted at point detectors in and near the emergent mouth of the duct are compared with 
similar results computed by the FASTER code and the OSR code using cylindrical geom- 
etry to describe the duct. The helical duct is given a very large pitch to simulate a 
cylindrical duct within the slabs considered. The comparison of number spectra is al- 
most exact, verifying the techniques developed herein. FORTRAN listings of computer 
subroutines that perform the numerical analysis are included; these subroutines are 
designed to be used in conjunction with FASTER, a high-efficiency Monte Carlo code. 


INTRODUCTION 

A useful reactor coolant flow passage geometry is that of a helicoid, or a helical 
tubular duct. This duct configuration has the dual advantage of reducing radiation stream- 
ing through a shield and permitting low coolant pressure drop. 

In order to track neutrons and gamma rays through such a complex geometric con- 
figuration, subroutines are written to be used with the FASTER code. Sample problems 

(ref. 1), previously run with cylindrical ducts with the FASTER and OSR codes, were 
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solved using a helical duct; this was done by assigning a pitch of 10 centimeters to a 
helicoid penetrating a 20- and a 100- centimeter-thick water slab. Thus, the important 
within-slab region of the helicoid is essentially cylindrical. 

Agreement was obtained for the neutron number spectrum at several point detectors 
in and near the emergent mouth of a duct having a point fission neutron source at the 

*ART Research Corporation, Los Angeles, California. 



center of one end of the duct. This agreement verified the helical duct geometry handling 
subroutines . 


ANALYSIS 


The equation of the helicoid is developed by using the fact that the intersection of the 
surface and a cutting plane normal to the centerline of the helicoid is circular. The pa- 
rametric equations of the duct centerline are 


x = R cos u> z 

(la) 

y = R sin o>z 

(lb) 


where w = 277/d, d is the pitch of the helical centerline, and R is the perpendicular 
distance from the axis of the cylinder around which the helicoid is wound to the helical 
centerline . 

A point (the origin of the double primed coordinate system is sketch (a)) is first fixed 
at an arbitrary azimuthal angle 9, Solving for the equivalent z-coordinate gives 
z Q = e/w. A coordinate system is centered at this point so that its y'-axis is tangent to 
the helical centerline and its x’-axis is along the cylindrical radius vector as in sketch (a): 



In the primed coordinate system, the intersection of the unwound helical duct with 
the cutting plane has the following equation: 

x’ 2 + z' 2 = a 2 (2) 

where a is the inner radius of the duct. Performing a rotation cp about the x'-axis, 
the appropriate transformation equations become 
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x’ = X M 


(3a) 

y' = y" cos (p + z" 

sin cp 

(3b) 

z* = -y” sin (p + z M 

cos <p 

(3c) 

sin (p = - 

1 

(3d) 

r 2 / . 21 1 / 2 / 

' 2 2\^/ 2 


|<f + (27rR)^J ( 

^1 + wTT J 


cos <p ^5: - 

wR 

(3e) 

r,2 / k21^/ 2 

/ o 2\*/ 2 


[cT + (2irR) j 

(l + coTr) 


x" 2 + (-y" sin <p + z" cos <p) 2 = a 2 

(4) 


The double primed coordinate system must also be translated and rotated as indi- 
cated in sketch (b) and equations (5). 

y" 

X" 

y 

(b) 



x” = x cos 9 + y sin 9 - R 

(5a) 

y ,s + -x sin 9 + y cos 9 

(5b) 

z" = z - z„ 
0 

(5c) 


so that 

fx cos 9 + y sin 9 - + f(-x sin 9 + y cos 9) sin q> + (z - z Q ) cos <p| 2 = a 2 (6) 
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To put everything in terms of the x, y , and z coordinates , the following relations 
are used: 



The resulting equation is 



( 7 ) 


(8a) 


(8b) 


It should be noted that the arctangent function subroutine will return only the principal 
value. Therefore, in the actual computations the quantity [tan"*(y/x) - wz] will be 
treated as modulo ti. 

Equation (8b) is generalized to the following form: 



( 9 ) 


where <p( r) is the helicoid surface function and has the following values: 
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<p(r) = 0 on the helicoid surface 
> 0 outside the surface 
< 0 within the surface 

andwhere (i,j,k) = (1, 2, 3), (2, 3, 1), or (3,1,2); = x,y,z; and x?,x°,x£ are 

coordinates of the origin of a translated coordinate system. 

Thus, the following data are normally supplied to describe the helical duct: 
i,x?,x°,x°,d,R,a. 

VALUES OF THE SURFACE FUNCTION 

The computer subroutine has been written to evaluate the surface function at the 
point r = r + s£2 where r Q , £2, and S are arguments supplied to the subroutine. Let 

r 0 = (x*,y*,z*) = (xj,x|,x|) 

£2 = (a,p,y) = (CpCgjCg) 

Then 

r = (x,y, z) = (x-pXgjXg) = (x| + CjSjX^ + C 2 S,x| + CgS) (10) 

The value <p(r) of the surface function is obtained by substituting 7 from equation (10) 
into equation (9). 

NORMAL VECTOR 

The components of the normal vector n are 

n = (n 1 ,n2,n 3 ) (11) 
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and are calculated as 



and 


l = 1,2,3 

The partial derivatives are calculated as 

= -2cuhv 
0x i 

A A 

».,(Y hvx k\ 

3x j \ ? ? 2 / 

8x k V ? ? 2 / 


where 


( x i- x i- 






h = 


R 2 

1 + oj 2 R 2 


( 12 ) 


(13a) 

(13b) 

(13c) 


(14a) 

(14b) 

(14c) 

(14d) 
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Intersections With Straight Line 


The intersection of a straight line with the helicoid is calculated by an iterative pro- 
cess. The assumption is made that the surface function is locally quadratic. Then the 
value of the function can be expressed by a truncated Taylor series as 


cp(r') = cp(r + SS2) = <p(r) + S 


3<?(r) + s 
as 


2 

2 



(15) 


where <p(r) is determined from equation (9), and the coefficient 0q?( r)/9S is calculated 
from 


3<p(r) = 

as 


0£(r) 


m 


3x. 


m 


using equation (13) to compute the partial derivatives. 
The final coefficient is given by 


as 2 



1=1 m=l 


c c a <pjd _ 

Z m 3x 3x 
l m 


where 


d 2 <?(r) _ 


2hw 


3X ? 


dV(r) _ 2 

.2 


3x 


3 


' m k . , h l™i k . *k 


? 3 ? 2 




f 


9 2 <p(r) = 2 

3x, 2 


~'2 "2 / „ ' ~ " 2 ' 

uxf x? /-2vx.x x. 

__I + _± + h L±+_i 


r 




,4 


C 4 


(16) 


(17) 


(18a) 


(18b) 


(18c) 
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(18d) 


d\(r) _ 2ha,x k 
3x^3Xj ^2 


d 2 (p(r) 

e Xj 3x k 


= 2 



(18e) 


(18f) 


The net result is a quadratic equation for the distance S to the surface since the 
surface function <p( r) must be zero at the intersection 


cp(r') = cp(Z) +sMI) + s! IjP.0 = Q 

as 2 as 2 


(19) 


Let 


Then 


A _ d 2 <?(r) 

as 2 

B _ 8<p(r) 
as 

C = 2<p(r) 


AS 2 + 2BS + C + 0 


(20a) 


(20b) 

(20c) 


(21a) 


or 


S = 


-B +6 


Vb 2 - AC 
A 


(21b) 


The sign of the square root is determined by the ambiguity index 5 (of eq. (21b)) 
with respect to the region in which the ray tracing is performed. The ambiguity index 
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is defined as 


-jfj 


( 22 ) 


where r g is any point in the region. 

Equation (21) is not used when | AC/B^| < 10"^. A one-term expansion is then used 
for the square root so that 


S 


£(,_ 6 B\ 1 _C5 

a\ M/ 2 N 


(23) 


This procedure eliminated convergence problems caused by underflow in the vicinity 
of the intersection. 

The distance S to the intersection is used to compute a new reference point r and 
the entire procedure is then repeated until the change in S is less than 10“^ of the value 
already obtained. In general, only two iterations are needed. 

The changes required to incorporate the helical duct in FASTER are shown in the 
appendix. The input instructions are modified to include the following on card 2-2: 

i , 0 , 16 , x° , y° , z° , L , R , a i th surface has helix axis parallel to x-axis 

i, 0, 17,y°, z°,x°, L,R, a i th surface has helix axis parallel to y-axis 

i,0,18,z°,x°,y°,L,R,a i*'* 1 surface has helix axis parallel to z-axis 


RESULTS AND DISCUSSION 

The duct analysis routine was verified by using it for the 20- and 100- centimeter- 
thick water slab cylindrical duct problem previously run as shown in figure 1 (ref. 1). 
Neutron number spectra essentially coincided with previously calculated spectra (05R 
and FASTER with cylindrical duct geometry) as seen in figure 2. The cylindrical duct 
was simulated with the helicoid geometry subroutines by making L, the pitch, equal to 
10 centimeters. Neutron number spectra at each detector point were calculated using 
1000 histories. Computer time was 20 minutes per detector point,, approximately double 
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that required for solution of the same problem using the cylindrical geometry routine , as 
expected because of the more complex geometric description of the cylindrical duct, 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, June 2, 1969, 

126-15-01-03-22. 
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APPENDIX -LISTING 

Changes to FASTER for Subroutine HELIX 


TO INCORPORATE SUBROUTINE HELIX IN FASTER* MAKE THE FOLLOWING CHANGES 
SUBROUTINE LOCDUM 
INSERT AFTER FAST0151 

IF(MAX.GT.g) CALL HELIX ( 1 *MAX~9* A ( 1 »K) »X»X*0«0*U(K) ) 

IF(MAX.GT,9> 60 TO 160 

SUBROUTINE GEOVIN 

REPLACE FAST0532 BY 

GO TO ( 232 * 233 * 235 * 260 , 255 , 282 ) * NGT 
INSERT AFTER FAST0579 
60 TO 290 
282 DO 284 J=l*6 
284 A(J*I) = AA(J) 

NTP(I) = NEX + 9 
INSERT AFTER FAST0616 

IF(MAX.GT,9> CALL HELIX ( 1 *MAX-o* A ( 1 *K ) * ADM, ADM, 0 ,0 ,FST) 

IF (MAX *6T ,9) GO TO 335 
INSERT AFTER FAST0619 
335 CONTINUE 

SUBROUTINE TRADUM 

INSERT AFTER FAST2228 

IF(NTP(K),LE,9) GO TO 198 
IF(M.GT.l) GO TO 192 
SI < 1 *K ) = STT 
SI <2*K ) = STT + SP 

CALL HELIX (3*KK*(NTP <K>-9) »A(1,K> *X»C*SI(1*K)*SR) 

IF<S8) 340*330*330 

192 CALL HELIX(1*NTP(K)-9,A(1*K)*X,C*STT*FST) 

IF < FLOAT (KK)*FST) 340,340*390 
198 CONTINUE 

SUBROUTINE NORDUM 

INSERT AFTER FAST2343 

IF(MAX,GT,9) CALL HELIX (2 , MAX-9* A f 1 * I) *X,X, 0 ,0 ,C ) 

IF(MAX.GT.9> GO TO 35 
C INSERT AFTER FAST2355 
35 CONTINUE 


FAST0151 

FAST0151 


FAST0532 

FAST0579 

FAST0579 

FAST0579 

FAST0579 

FAST061.6 

FAST0616 

FAST0619 


FAST2228 

FAST2228 

FAST2228 

FA5T2228 

FAST2228 

FAST2228 

FA5T2228 

FAST2228 

FAST2228 


FAST2343 

FAST2343 

FAST2355 
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Subroutine HELIX Listing 


C***********************************************************************hELIXOoO 


C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


(I.J»K)=(1,2,3)=(X.Y,Z) 
tI.J»K)a(2,3,l)={Y»Z.X) 
(I.J.K)=(3.1 »2)=(Z»X.Y) 


/JARG/ = 1* HELIX PARALLEL TO X-AXIS 

/JARG/ = ?. HELIX PARALLEL TO Y-AXTS 

/JARG/ = 3. HELIX PARALLEL TO Z-AXIS 

SIGN(JARG) = AMBIGUITY INDEX 

VARG(l) = REFERENCE VALUE OF X(I) AT WHICH THETA = 0,0 
VARG ( 2 ) a TRANSLATION OF HELIX AXIS IN X(J) DIRECTION 
VARG ( 3 ) = TRANSLATION OF HELIX AXIS IN X(K) DIRECTION 
VARGC4) = PERIOD OF HELIX ALONG X Cl) AXIS 

VARG (5) = RADIUS OF GYRATION OF HELIX PERPENDICULAR TO X(I) AXIS 

VARG (6) = RADIUS OF HELIX CROSS SECTION 

WARG(l) S X COORDINATE OF ORIGIN OF RAY 

WARG(2> = Y COORDINATE OF ORIGIN OF RAY 

WARG<3> S Z COORDINATE bF ORIGIN OF RAY 

XARG ( 1 ) = DIRECTION COSINE OF RAY WITH RESPECT TO X AXIS 

XARG ( 2 ) = DIRECTION COSINE OF RAY WITH RESPECT TO Y AXIS 

XARG (3) = DIRECTION COSINE OF RAY WITH RESPECT TO Z AXIS 

YARG(l) = DISTANCE ALREADY TRAVERSED ALONG RAY 

YAR6 (?) 


MAXIMUM DISTANCE ALONG RAY 
I ARG r 1, COMPUTE VALUE OF HELIX 
IAR6 = 2. COMPUTE NORMAL VECTOR 
I ARG a 3 » COMPUTE DISTANCE TO HELIX 


PETURNED IN ZARGC1) 

RETURNED IN 7 ARG ( 1 ).,,,. 7ARp { 3) 
RETURNED IN ZAPG(l) 


HELIXOOO 

HELIXOOO 

HELIXOOO 

HELIXOOO 

HELIXOOO 

HELIXOOO 

HELIXOOO 

HELIXOOO 

HELIXOOO 

HELIXOOO 

HELIXOOO 

HELIXOOO 

HELIXOOO 

HELIXOOO 

HELIXOOO 

HELIXOOO 

HELIXOOO 

HELIXOOO 

HELIXOOO 

HELIXOOO 

HELIXOOO 


C**********4^4^*********************************************************HELTX000 


SUBROUTINE HELIX ( I ARG. JARG .VARG .WARG .XARG. YARP. ZARG) 
DIMENSION VARG ( 1 ) » WARG ( II » XARG ( 1 ) . YARD ( 1 ) . ZARG ( 1 ) 
DIMENSION XZ(3) »X(3) ,C(3> .CN(3) 


hElixoio 

HELIX020 
HEL 1X030 


DATA 

KOUNT.PI.TPI/0. 3. 1415927.6, 2*31853/ 

HELTX040 

DO 100 

1=1.3 

HELIX050 

cm 

= XARG ( I ) 

HELIX060 

xzm 

= WARG ( I ) ♦ Cm*YARG(l> 

HELIX070 

100 X ( I ) 

= XZ(I) 

HELIXOOO 

NGT 

= IARG 

HFLTY090 

I 

= IABS(JARG) 

helixioo 

J 

= I + 1 - 3* ( T/3) 

HFLIX110 

K 

= J + 1 - 3*(J/3> 

HELTX120 

TPL 

= TPI/VARG(4) 

HELIX130 

GAM 

= 1 »0/(TPL**2 + 1,0/VAPG(5)**2) 

HELIX140 

110 DXI 

= X ( I ) - VARG ( 1 ) 

HFLIX150 

DXJ 

= X(J) - VARG ( 2 ) 

HELIX 160 

DXK 

= X (K ) - VARG(3> 

HEL IX 170 

USM 

= SORT (DXJ**2 + DXK**2) 

HELIX180 

THE 

= 0.0 

HELIX185 

IF (USM, 

GT.0,0) THE = ATAN2(DXK,DXJ) - TPL*DXI 

HELIX190 

VSM 

= THE - TPI*FLOAT ( IFIX ( (THE + PD/TPD) 

HELTX200 

UBG 

= USM - VARG(5) 

HELIX210 

PHI 

= UBG**2 + 6AM*VSM**2 - VARG(6>**2 

HELIX220 

IF(NGT, 

GT, 1 ) GO TO 120 

HELIX230 

ZARG(l) 

= PHI 

HEL 1X240 

GO TO 990 

HELIX250 

120 CST 

a DXJ/USM 

HEL 1X260 

SNT 

= DXK/USM 

HEL 1X270 

GNU) 

= -TPL*GAM*VSM 

HELIX200 

CN(J> 

a UBG*CST - GAM*VSM*SNT/USM 

HELIX290 

Cm ( k ) 

= UBG*SNT + GAM*VSM*CST/USM 

HEL 1X300 

OPH 

= 0,0 

HELIX310 

IF (NGT, 

GT,2) GO TO 150 

HELIX320 

DO 130 

1=1.3 

HELIX330 

130 OPH 

= DPH + CN(I)**2 

HELIX340 

DPH 

= SORT (DPH) 

HELIX350 

DO 140 

1=1.3 

HELIX360 

140 ZARG(I) 

= CN ( I ) /DPH 

HEL 1X370 

GO TO 990 

HELIX30O 
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150 DO 160 
160 DPH 
CN(I) 
CN(J) 
CN(K> 
DDP 

DO 170 


L=l,3 

= DPH + C(L)*CN(L) 
s GAM*TPL**2 

= URG*SNT**2/USM+CST**24-GAM*SNT*(2.0*V5M*CST+SNT)/USM**2 
= UB6*CST**2/USM+SNT**2-GAM*C$T*(2.0*VSM*SNT.-eST)/USM**2 
= 0.0 
L=l»3 


170 DDP 
CN ( I ) 
CN(J) 

1 

CN(K) 
DO ISO 
M 

ISO DDP 

IF(NGT 
TOT 
KK 
SGN 
STB 
190 H 


= DDP + CN(L)*C(L)**2 
= TPL*GAM*SNT/USM 
= SNT*CST* ( 1.0 - UBG/USM) 

+ GAM* ( VSM* ( 5NT**2 - CST**2) 
= -TPL*GAM*CST/USM 
L=l,3 

= L + 1 - 3*(L/3) 
s DDP + 2 , 0 *CN ( L ) *C ( L ) *C ( M ) 
»GT,3) GO TO 190 
= 0.0 
= I/JARG 
= FLOAT (KK) 

= 2 »0*YARG (2) 

= 0.0 


SNT*CST)/USM**2 


NFH = 0 

DLS = 0.0 

IF (PHT .NE.O.O) GO TO 200 
IF (DDP) 210*250.210 
200 IF(DDP. NE.O.O) GO TO 210 
DLS = -0,5*PHI/DPH 

GO TO 250 

210 IF(DPH.EG.O.O) GO TO 220 

IF(ARS(DDP*PHI/DPH**2) .GT.l.OF-4) GO TO 220 
ADP = ABS(DPH) 

DLS = ADP* (SGN - DPH/ADP)/DOP - 0 ,5*SGN*PHT/ADP 

GO TO 250 

220 DLS = -DPH/DDP 

H = DPH**2 - DDP*PHT 

IF(H) 230.250.240 
230 NFH = 1 

GO TO 250 

240 DLS = DLS + SGN*SORT (H) /DDP 

250 I F ( (DLS + TOT). LF. 0,0) GO TO 290 

TOT = TOT + DLS 

IF ( (NFH.GT.O) .AND, (NGT.GT.3) ) GO TO 290 
IF(NFH.GT.O) GO TO 260 
IF ( ABS (DLS/TOT) «LF.1»0E”5) GO TO 280 
IF (NGT ,GT,4 ) GO TO 280 
260 DO 270 L=l*3 
270 X (L) * XZ(L) + TOT*C (L) 

NGT s NGT + 1 

GO TO 110 

280 STB = TOT 

IF(DLS.LE.O.O) STB = STB + AMAX1 (-OLS,1.0E-5*TOT) 
290 ZARG(l) = STB 


KOUNT = KOUNT + 1 

IF (K0UNT.LE.1D0 ) WRITE (6. 2000 )NGf .I.J.K.KK.NEH.X.C, DDP. DPH, PHI, 
1T0T , DLS, STB 
990 RFTURN 

2000 FORMAT ( 1X»6I6,1P6E12.4/37X,6E12»4) 
end 


HELIX390 
HFLTX400 
HFLIX410 
HEL 1X420 
HELIX430 
HFLIX440 
HELTX450 
HELIX460 
HELIX470 
HELTX480 
HELIX485 
HEL 1X490 
HELTX500 
HELTX510 
HELIX520 
HELIX530 
HELTX540 
HFLIX550 
HPLTX560 
HELIX570 
HELIX580 
HELIX590 
HELIX600 
HELIX610 
HELIX620 
HEL 1X630 
HFLIX640 
HEL 1X650 
HELTX660 
HEL 1X670 
HELTX680 
HELTX690 
HELIX700 
HEL 1X7 10 
HELIX720 
HELTX730 
HELIX740 
HELIX750 
HEL 1X760 
HELTX770 
HELIX780 
HEL 1X790 
HELTX800 
HELIX810 
HELIX820 
HEL1X830 
HEL 1X840 
HELIX850 
HEL 1X860 
HELIX870 
HELIX880 
HELIX890 
HEL 1X900 
HELIX910 
HELTX920 
HEL 1X930 
HEL 1X940 
HELIX950 
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(a) Detector 1; slab thickness, 20 centimeters. (b) Detector 2; slab thickness, 20 centimeters. (c) Detector 4 ; slab thickness, 20 centimeters. 
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